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Transient dynamics is pervasive in the human brain and poses chahenging problems both in mathematical 
tractability and clinical observability. We investigate statistical properties of transient cortical wave patterns 
with characteristic forms (shape, size, duration) in a canonical reaction-diffusion model with mean field 
inhibition. The patterns are formed by a ghost near a saddle-node bifurcation in which a stable traveling 
wave (node) collides with its critical nucleation mass (saddle). Similar patterns have been observed with 
fMRI in migraine. Our results support the controversial idea that waves of cortical spreading depression (SD) 
have a causal relationship with the headache phase in migraine and therefore occur not only in migraine with 
aura (MA) but also in migraine without aura (MO), i.e., in the two major migraine subforms. We suggest 
a congruence between the prevalence of MO and MA with the statistical properties of the traveling waves' 
forms, according to which (i) activation of nociceptive mechanisms relevant for headache is dependent upon 
a sufficiently large instantaneous affected cortical area anti-correlated to both SD duration and total affected 
cortical area such that headache would be less severe in MA than in MO (ii) the incidence of MA is reffected 
in the distance to the saddle-node bifurcation, and (iii) the contested notion of MO attacks with silent aura 
is resolved. We brieffy discuss model-based control and means by which neuromodulation techniques may 
affect pathways of pain formation. 



features are central and are by no means exclusive to bi- 
ological membranes but shared by all excitable elements. 
Firstly, the inevitable threshold in any all-or-none behav- 
ior requires nonlinear dynamics. Secondly, the transient 
response of the system to a super-threshold stimulation 
eventually has to lead back to a globally stable steady 
state after some large phase space excursion. This in- 
dicates global dynamics, that is, dynamics involving not 
only fixed points and their local bifurcations but larger 
invariant sets, for instance periodic orbits that collide 
with fixed points. An excitable element is in some sense 
the washed-up brother of the relaxation oscillator: when 
the threshold vanishes, a single excitable element usually 
becomes a simpler behaved — and much longer known — 
relaxation oscillatoi^. 

In this study, we propose a model for wave patterns 
with a characteristic shape, size, and duration. These 
waves are transient responses to confined, spatially struc- 
tured perturbations of the homogeneous steady state. 
The homogeneous steady state is globally stable in our 
proposed reaction-diffusion model because we also intro- 
duce an effective inhibitory mean field feedback control. 
This leads to a new type of local excitability in a spa- 
tially extended medium involving disappearing traveling 
wave solutions as larger invariant sets. Both, the model 
and the initial conditions are motivated by thepatho- 
physiology of migraine and clinical observations^^^ and 
the results are applied to some currently controversial 
topicP^l. 

We will briefiy introduce concepts of excitable elements 
and excitable media in two-variable reaction-diffusion 
systems and also the idea of an additional long-range in- 
hibitory feedback that is studied in various other systems 

outside the neurosciences and also in neural field models. 
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AUTHOR SUMMARY 

The key to the genesis of migraine with aura is a 
traveling wave phenomenon called cortical spreading de- 
pression (SD). Migraine is characterized by recurrent 
episodes, the aura phase usually lasts only about 30 min. 
Thus, SD is a transient state. During its course, SD 
massively perturbs the brain's ion homoeostasis. We re- 
solve the puzzling problem why SD does not engulf all 
of the densely packed excitable neurons by suggesting 
a well-established pattern formation mechanism of long- 
range inhibitory feedback. Furthermore, we use cortical 
feature maps to create plausible initial conditions as per- 
turbations of the homogeneous state. The statistics of 
occurrences of the different classes has the potential to 
reproduce epidemiological statistics of different diagnos- 
tic forms of migraine such as migraine with or without 
aura and provides simple answers to some very controver- 
sially discussed current questions in migraine research. 



INTRODUCTION 

The undoubtedly most fundamental example of tran- 
sient dynamics is the dynamical phenomenon of excitabil- 
ity, that is, all-or-none behavior. Shortly after transient 
response properties of excitable membranes were clas- 
sified into two classed, it was also explained in a de- 
tailed mathematical model how excitability emerges from 
electrophysiological properties of such membranes in the 
ground-breaking work by Hodgkin and Huxle}/^. Two 
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graine as a dynamical disease is more elaborated in the 
discussion. 

The original conductance-based membrane model from 
Hodgkin and Huxley, and the more refined versions to 
date, contain many variables, but fortunately this is not 
essential for excitable elements. In fact, it turned out 
that the two classes of excitability are actually amenable 
to direct analysis in a two-dimensional phase plane by 
identifying fast and slow processes in the conductance- 
based model and grouping these into dynamics of just 
two lump variables^^ Using such a geometrical ap- 
proach and partly analytical theory, the original empiri- 
cal classification of excitability was further pursued with 
bifurcation analysiJ^, explaining class I by identifying 
its threshold as a stable manifold of a saddle point on 
an invariant cycle and the one of class II as a trajectory 
from which nearby trajectories diverge sharply (called a 
canard trajectory). Extensions to these principal mecha- 
nisms involve codimension 2 bifurcations and lead also to 
bursting in three- variable models, which have been inves- 
tigated in great detail^^. However, the two- variable mod- 
els of a fast activator and slow inhibitor and their phase 
portraits of class I and II became qualitative prototypes 
for excitable systems in various biologicaf^, chemicaP^ 
and physical contexts^^. 

Distinct from these excitable systems are those that are 
spatially extended systems, called excitable media. (We 
use the word "system" as a general term, "medium" only 
for spatially extended systems, and "element" for point- 
like systems.) Already the original work by Hodgkin and 
Huxle}El described extended tube-like membranes (ax- 
ons) and introduced the cable equation as a parabolic 
partial differential equation, which is in the same class as 
the diffusion equation. Even in reaction-diffusion media 
with infinite-dimensional phase space, we can again apply 
geometrical approaches, simply because excitable media 
are not defined by — in contrast to excitable elements — 
transient dynamics but traveling wave solutions. The 
quiescent state is the non-excited homogeneous steady 
state. And like the quiescent state, excited states of a 
medium are usually stationary states in some appropriate 
comoving frame with ^ = x—ct. Furthermore, the thresh- 
old is related to an unstable stationary state, the critical 
nucleation solution, usually in another comoving frame 
including c = 0. The existence of the nucleation solution 
is a simple consequence of multistability, see Fig.[l]A, but 
note that even monostable excitable elements have simi- 
lar unstable stationary states in class I. 

Central to our approach is a distinct subexcitable 
medium in which localized traveling waves occur only 
transiently. Reaction-diffusion waves would engulf all 
of the medium, if formed in a two-variable system with 
only one activator and one inhibitor with the system's 
parameters in the appropriate regime. In contrast, lo- 
calized traveling waves indicate a demand-controlled ex- 
citability. Similar ideas to obtain localized traveling, 
though not transient, waves have been introduced in 
various contexts, for instance an integral negative feed- 




FIG. 1. Schematic sketch of the phase space of (A) the 
uncontrolled system, (B) the system with mean field control 
adjusted such that the nucleation solution is stabilized and 
(C) the system with mean field control and control parame- 
ters such that the ghost of the saddle-node bifurcation is still 
influencing the dynamics. 



back or a third, fast diffusing inhibitory component for 
moving spots in semiconductor materials, gas discharge 
phenomena, and chemical system^^^UHl^ Furthermore, in 
neural field models^^, localized two-dimensional bumps 
are studied^^^ in integrodifferential equations (with- 
out diffusion) in the context, for example, of memory 
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formatiorP^. Localized structures have also been dis- 
cussed in the context of cortical spreading depression 
(SD) in migraine before, in particular a model wit h nar- 
rowly tuned parameters that shows transient 
and a model with mean field feedback control that allows 
for localized waves^^. But it is for the first time now that 
a model is presented in which wave phenomena occur that 
are both localized and transient, so that a variety of new 
questio ns that a re controversially discussed in migraine 
can be adressed. 
Migraine is characterized by recurrent episodes of head 
pain, often throbbing and unilateral. In migraine without 
aura (MO), attacks are usually associated with nausea, 
vomiting, or sensitivity to light, sound, or movementP^. 
Migraine with aura (MA) involve, in addition but also 
rarely exclusively, neurologic symptoms (aura) that are 
associated with waves of cortical SE0^. SD is a reaction- 
diffusion process, although, clearly, the originally pro- 
posed mechanism of a simple one-variable model for 
SD front propagation triggered by diffusion of elevated 
extracellular potassium with bistable kinetics, to date 
known as Hodgkin-Huxley-Grafstein model^^, does not 
capture the complex chain of involved reaction j^^^^^l 
However, migraine aura symptoms manifest themselves 
on a macroscopic scale over several minutes up to one 
hour and extend over several centimeters when mapped 
onto the corresponding cortical areas^^^, see Fig. 2 In 
this study, we are interested in these clinically relevant 
properties of migraine with aura. To this end, we ex- 
ploit the concept of nucleation, growth, and subsequent 
shrinking of SD in a canonical reaction-diffusion model of 
activator-inhibitor type with mean field feedback control. 



RESULTS 

We first present a model constructed by adding a 
mean field inhibitory feedback to a well known reaction- 
diffusion system, then we present the statistical proper- 
ties of the transient behavior this model is exhibiting. 



Before infinity and beyond by mean field inhibitory 
feedback control 

The diversity of the behavior of traveling waves in 
two spatial dimensions was studied in canonical models 
(see Discussion) depending on the two generic param- 
eters /3 and £ in Eqs. ([T])-([2|), which determine the pa- 
rameter plane of excitabilit}'^''. Eqs. ([l])-([2| determine an 
excitable medium without feedback control. In those me- 
dia, patterns of discontinuous (open ends) spiral-shaped 
waves are used to probe excitability and these patterns 
are closely related to the discontinuous, localized tran- 
sient waves we propose in our model. 

We make use of the fact that at a low critical excitabil- 
ity, called the rotor boundary dRoo^ spiral waves do not 
curl-in anymore but become half plane waves^^^^. Be- 



yond the rotor boundary lies the subexcitable regime in 
which discontinuous waves start to retract at their open 
ends and any discontinuous wave is transient and will 
eventually disappear. The border dRoo marks a saddle- 
node bifurcation at which discontinuous waves collide 
with their corresponding nucleation solution. This leads 
to the key idea of our model. A linear mean field feedback 
control moves this saddle-node bifurcation towards dis- 
tinct localized wave segments with a characteristic form 
(shape, size) and behind this bifurcation these waves be- 
come transient objects (see Fig. [T]). 

Before we introduce the effect of mean field feedback 
control, we have to consider the behavior of continuous 
waves (closed waves fronts without open ends) when the 
excitability is decreased. This will be important if we 
want to understand the fate of any solution, discontinu- 
ous or not, under mean field feedback control. Unbroken 
plane waves propagate persistently even if the parame- 
ters are chosen in the subexcitable regime until the prop- 
agation boundary dP is reached. At this border, the 
medium's excitability becomes too weak for continuous 
plane waves to propagate persistently. The border dP 
in parameter space indicates again a saddle-node bifur- 
cation at which a planar traveling wave solution collides 
with its corresponding nucleation solution. Note, that 
the planar wave is essentially a pulse solution in ID and 
the nucleation solution in ID is called the slow wav^^. 

In Fig. |3]A., both the rotor boundary dRoo and the 
propagation boundary dP are shown in a bifurcation di- 
agram for the excitable medium described by Eqs. ([T])-([2|. 
We chose j3 as the bifurcation parameter and follow (see 
Methods) the branch of the unstable nucleation solution 
(NS) whose stable manifold separates the basins of at- 
traction of the homogeneous state and a spiral wave (with 
two counter-rotating open ends). The unstable manifold 
of NS consists of the two heteroclinic connections, one to 
the stable homogeneous state and the other to the trav- 
eling wave solution (see Fig.jl]). The order parameter on 
the ordinate in Fig. |3] is the surface area S inside the 
isoclines at = of the traveling wave solutions, see 
Eq. (|3| . We call S the wave size. 

The mean field control that we introduce by Eqs. ([3|- 
Q establishes a linear feedback signal of the wave size 
S to the threshold j3. With this linear relation, we in- 
troduce two new parameters, the coupling constant K 
and /3o, the threshold parameter for the medium without 
an excited state {S = 0). Note that the parameter /3o 
can be also seen as the sum of two threshold values, the 
former (3 in Eq. ^ and an offset coming from the new 
control scheme. While the introduction of the control 
introduces two new parameters /3o and i^, at the same 
time 13 becomes dependent upon the control, so that we 
have a total of three parameters, in this study, we kept 
K fixed dX K = 0.003 and varied /3o, with a particular 
focus on the statistics for /3o G [1.32,1.33,1.34]. 

We chose /3o as the new bifurcation parameter in the 
bifurcation diagram for the full reaction-diffusion model 
with mean field coupling described by Eqs. Q-Q, see 
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FIG. 2. Source localization of the magnetic resonance (MR) data signal of SD (from'^. Color code: time from onset, locations 
showing the first MR signals of SD are coded in red, later times are coded by green and blue (see color scale to the right). 
Signals from the first 975 seconds were not recorded because the migraine attack was triggered outside the MR imaging facility. 
(A) The data on folded right posterior pole hemispheric cortex; (B) the same data on inflated cortical surface; (C and D) the 
same data shown on the entire hemisphere from posterior- medial view (oblique forward facing), folded and inflated, respectively. 
As described in the original studjA, MR data were not acquired from the extreme posterior tip of the occipital pole (rearmost 
portion). (E) A fully flattened view of the cortical surface. The aura-related changes are localized wave segments. Note that in 
the flattened cortex was cut along the steep sulcus calcarine to avoid large area distortions induced by the flattening process. 
The colored border to the left is the cut edge that should be considered being connected such that the color match up as seen 
in (B and D). (Copyright permission from authors granted, from PNAS is requested.) 



Fig. |3^. This diagram is a sheared version of the one 
without mean field coupling in Fig. [3]A. While it is a 
trivial fact, that the linear relation in Eq. Q describes 
an affine shear of the axises (/3, 5) of bifurcation dia- 
gram in A to the axises (Po^S) in B, the fact that the 
branch of the NS solution can be mapped this way is not. 
Firstly, this relies on the way we introduce the feedback 
term. It just adds a constant value to the old bifurca- 
tion parameter /3, if the solution under consideration is 
stationary. Therefore, any stationary solution must ex- 
ist in both diagrams being just sheared branches. The 
same holds true for traveling wave solutions that are sta- 
tionary in some appropriate comoving frame £^ = x — ct 
with speed c. However, not much can be said about the 
stability of such solutions, when we introduce the mean 
field feedback term. 

The branch of the formerly unstable nucleation solu- 
tion NS (Fig. ^) folds in Fig. such that two solu- 
tions coincide for a given value of until they collide 
and annihilate each other at a finite value of 5 ~ 5.5 for 



K = 0.003. For the fixed value of K = 0.003, the upper 
branch is a stable traveling wave solution in the shape 
of a wave segment, while the lower branch is the corre- 
sponding nucleation solution of this wave segments, as 
schematically shown in Fig. [l^. The fact that the upper 
branch is stable was confirmed by numerical simulations. 
Larger that is, a less steep control line in Fig. [3|^ 
can be seen as a "harder" control, because a small given 
change in S leads to larger variations in the effective pa- 
rameter /3. As a consequence, it is difficult to continue by 
means of this control the lower part of the branch corre- 
sponding to small traveling wave segments in numerical 
simulations. 

The choice of a parameter regime for this model that 
shows transient localized waves and is globally stable 
with the homogeneous state as the only attractor is now 
straightforward. Transient localized waves occur due to 
a bottleneck — or ghost behavior — after the saddle-node 
bifurcation. 
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FIG. 3. (left) The S - {5 plane with nucleation solution NS , propagation boundary ^P, rotor boundary and control 
lines for the used values of /3o. (right) The S - /So plane with the same quantities. 



Statistical properties 




FIG. 4. An example of a transient solution. Initial condi- 
tions, i.e. activator concentration u at Os (upper left), snap- 
shot of activator-concentration u after 30 s (upper right), after 
180 s (lower right). Time of passing through threshold value 
1^0 = from below the first time, i.e. passing of the wave 
front (lower left). 



To examine the typical transient patterns that the sys- 
tem generates, we want to know how the system responds 
to arbitrary initial conditions with the noticeable con- 
straint that the system should initially be in the homoge- 
neous steady state almost everywhere and the arbitrary 
perturbations from it are localized. As it is not possi- 
ble to formulate an analytical solution to the equations 
for arbitrary initial conditions, the idea is to simulate 
the dynamics for (many) different initial conditions. Ide- 
ally, these initial conditions should be "equally spaced" in 
phase space, in order to obtain relevant statistics about 
the different evolution possibilities. The problem that 
arises at this point is that an initial condition of this sys- 
tem is not only living in an infinite dimensional space 
(an initial condition would be given by two -functions 
R) but because of the nonlinearity of the equations, 
the set of all solutions is not even a vector space. To our 
knowledge, there is no helpful mathematical structure 
that could guide us in choosing our initial conditions. To 
attack this problem, we take a set of patterns, which are 
parameterized by a finite number of parameters and scan 
through these parameter. For details on the patterns see 
Methods. 

Of course, in characterizing the solutions, the same 
problem arises and appropriate characteristic parameters 
for the solutions have to be defined. 

To explain the three parameters we have chosen for the 
solutions and why they suit this problem, it is helpful to 
have a look at the lower left part of Fig. |4j in which an 
example solution is displayed. The first parameter we 
choose is the maximal area in which such a solution has 
activator concentration over a certain threshold level at 
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Distribution of solutions for the control close to 



FIG. 6. Distribution of solutions for the control line at 
an intermediate distance to ^r. 



one instant of time, termed maximal instantaneous area 
(MIA). The threshold level is taken toheu = 0, although 
this is the same threshold as used to define 5*, this is 
rather convenience than necessity. The second param- 
eter is the total area that has experienced an activator 
concentration above this level at some time during the 
course of the solution, termed total affected area (TAA). 
The third parameter is the time, during which the area 
of activator concentration above threshold is non-zero, 
termed the excitation duration (ED). Of course, the ex- 
act value of all these parameters for one single solution 
depends on the choice of threshold. For once, the thresh- 
old value has to be chosen such that after the activator 
concentration has fallen below it, no secondary excitation 
will be generated. 

The example solution depicted in Fig. |4] is a compar- 
atively long lived solution. It starts out very symmetri- 
cally (circularly) shaped, at one instant of time it breaks 
open into a discontinuous wave and a shape of the front 
develops, which is similar to that of a particle-like wave 
but because of the chosen control parameters, it shrinks 
in time and vanishes in the end. Because at the point 
when the circle breaks open, a comparatively large area 
is affected, it takes some time until it vanishes and the re- 
sulting TAA is relatively large. So this example solution 
has large ED, large MIA and large TAA. If the circle had 
not broken open at all, the control would have made the 
threshold value very large and the solution would have 
collapsed very quickly because of the propagation bound- 
ary 9P, such that the ED and the TAA would have been 
short, whereas the MIA would have been large. Other 



prototypical courses of solutions take place for instant, 
when the initial conditions affect the activator over a 
larger area but only in the middle of the area, the value 
is high enough to start a solution. In the surrounding 
area, the activator level is not high enough for that but 
the increased activator concentration leads to a rise in 
inhibitor concentration until the time the front reaches 
those parts and as a consequence, the solution vanishes 
early, having small ED, small TAA and small MIA as a 
consequence. 

We did the simulation for three different adjustments 
of the control force, successively going farther and farther 
away from the bifurcation point. Each of these simula- 
tions were started using 8000 initial conditions gener- 
ated in a manner that is described in Methods. Each 
of these initial conditions resulted in a solution that was 
classified according to the three parameters mentioned 
above. The solutions that did not result in any exci- 
tation at all (ED=MIA=TAA=0) were discarded. The 
density plots according to the classification parameters 
are shown in Figs. [5j |6j jT] First of all, though the dis- 
tribution of the solutions varies significantly, the number 
of solutions that represent an excitation hardly varies at 
all (4171 for small, 4183 for intermediate, 4182 for large 
distance from the bifurcation point), the symmetric dif- 
ference between the sets of initial conditions that lead to 
an excitation contains between 5 and 19 solutions. From 
this we can also deduce that the set of initial conditions 
that lead to an excitation does not significantly depend 
on the choice of control parameters. 

When looking at Fig. |5] one notices a clustering of 
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FIG. 7. Distribution of solutions for the control far away 
from ^R. 



the solutions in certain regions of the classification pa- 
rameters. In the section that depicts the TAA against 
the MIA, we notice three coarse clusters. Cluster I, the 
largest with high MIA and comparatively low TAA; clus- 
ter II, one that is less populated with low MIA and low 
TAA; and cluster III, one that is very sparsely populated 
with intermediate MIA and high TAA. The boundaries 
between these clusters are not very sharp. One could 
think that a solution that affects an overall large area 
(high TAA) will also affect a large area at one instant of 
time (high MIA). From looking at the mentioned clus- 
ters, one sees that this is not the case, the solutions with 
the highest MIA have all comparatively low TAA (clus- 
ters I and II) and the ones that have a high TAA only 
achieve an intermediate MIA (culster III). 

A partial explanation for this can be read off from the 
depiction of TAA against ED. All solutions that have a 
large TAA are also solutions that have a large ED, i.e., 
cluster III is distinct also in this plane. More than that, 
the dependence seems to be almost linearly. This is rem- 
iniscent of the localized particle-like wave solutions. For 
these, the area that is affected grows linearly in time be- 
cause the area that these solutions occupy at one instant 
of time is constant. 

The two clusters I and II that we observed merge to 
one in this plane of projection because they differ only 
very little in ED. This can also be noted, when comparing 
the planes MIA vs. TAA and MIA vs. ED, also here the 
cluster III with high TAA translates to a cluster with high 
ED and the cluster I with high MIA and comparatively 
low TAA moves closer to cluster II with both low ED and 



MIA. 

When varying the /3o parameter of the control force, 
the distribution of solutions in MIA-TAA-ED-space 
changes drastically. Upon raising the /3o parameter from 
/3o = 1.32 over f3o = 1.33 to /3o = 1.34, the system is 
put more and more into the subexcitable regime and the 
solutions are less and less affected by the ghost behavior 
(saddle- node bifurcation), see Fig. [sj This is noticeable 
by observing that the cluster with high TAA / high ED 
becomes less pronounced and vanishes almost completely 
for Po — 1.34. This can be understood as an interplay 
between the mean value of MIA in cluster III at about 
25 and S at the propagation boundary (at 9P, S ^ 24, 
S ^ 20.75, and S ^ 17.5 for f3o = 1.32, f3o = 1.33, and 
/3o = 1.34, respectively). For the control line farthest 
away from the saddle-node bifurcation (/3o = 1.34), dP is 
below even the smallest values of MIA in cluster III. Note 
that the value of S at the ghost is about 6, well below the 
propagation boundary. Also the other two clusters merge 
though there still exist solutions with high and with low 
MIA, but the transition is much more fuzzy than it was 
before. 

In Figs. [5) [6] and [T) we have included a little 'bestiary' 
to illustrate the typical courses of solutions in the respec- 
tive clusters and their change upon varying the param- 
eter /3o, the initial conditions for solutions 1-4 in these 
figures are always the same. From this arbitrarily chosen 
selection, we see that the MIA of each solution hardly 
changes between the Po values. Whereas the change of 
TAA and ED always go hand in hand and — depending 
on the cluster — can be up to fourfold for the chosen range 
of /3o. 

One could argue that the formation of clusters is an 
artefact of the choice of initial conditions. There is no 
simple answer to this. As mentioned, it is not possible to 
examine the complete set of initial conditions. Neither 
does this set carry a helpful structure which would al- 
low a sensible 'equidistant' sampling. This is the reason 
why we made the mentioned choice of initial conditions. 
For testing purposes, we also tried different schemes for 
the generation of initial conditions and found the same 
distribution of clusters qualitatively. 

In Fig. [8) we have plotted the cumulative distribution 
functions for the three classification parameters and the 
three choices of mean field control. From this picture we 
see, that the distribution of the MIA is only hardly influ- 
enced by the choice of control. This is very different for 
TAA and ED. For the TAA for example there are values 
(around 75), where for one choice of control the majority 
of solutions is below and for another choice the major- 
ity is above. For example, the fraction of values below 
TAA=80 is 0.995 for Po = 1-34 and 0.216 for /3o = 1.32. 
Also, we see that the cumulative distribution function for 
the TAA converges to 1 much slower, the closer the con- 
trol is to the saddle-node bifurcation. This means, that 
more solutions with high TAA exist for these choices of 
control. 
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FIG. 8. Cumulative distribution functions for the different classification parameters. 



DISCUSSION 

In this section we discuss three subjects related to 
the intended application to migraine pathophysiology. 
Firstly, the possible congruence between the prevalence 
of migraine subforms with the statistical properties of the 
wave patterns we observed; secondly, the possible phys- 
iological origin of the inhibitory feedback control; and 
thirdly, novel therapeutic approaches. We start with a 
discussion of the approach of using a canonical model. 



Canonical model and generic parameters for weakly 
excitable media 



More realistic models of SD are given by a 
conductance- ion- based models of SlW^ with up to 29 
dynamics variables. The fast inhibitory feedback that 
we suggest can be modeled in addition by neural field 
models^^. Such large-scale models of brain structure, in- 
cluding lateral connectivity, are available but still require 
enormous computer capabilities. We will argue here, in 
which sense our model is canonical for the problem we 
attack. 

Generally speaking, an excitable medium is a spatially 
extended system with a stable homogeneous steady state 
being the quiescent state and one or many excited states 
that develop after a sufficient perturbation from the qui- 
escent state (Fig. [l}^). The excited states are traveling 
wave solutions that propagate with a stable profile of 
permanent shape (possibly with some temporal modula- 
tion, such as breathing or meandering). To study generic 
features of an excitable medium, the simulations are of- 
ten carried out in the reaction-diffusion system given 
by Eqs. ([T])-([2|, the popular FitzHugh-Nagumo kinetics. 
Originally, the FitzHugh-Nagumo kinetics were a cari- 
cature of t he ele ctrophysiological properties of excitable 
membrane d^^ l ^^ l, but these equations with D = became 
a canonical model of local excitability of type II (based 



on Hopf bifurcation, either supercritical with subsequent 
extremely fast transition to a large amplitude limit cycle, 
named canard explosion^, or subcritical^^), and also, for 
I) ^ 0, of spatial excitabilit}^^, sometimes including dif- 
fusion in the second inhibitory species, which we do not 
consider here. Because we are considering transient be- 
havior originating from a high threshold regime (towards 
weak excitability), the classification of local excitability 
in type I and II (based on the transition at vanishing 
threshold, i.e., into the oscillatory regime) is not rele- 
vant. Furthermore, it is not clear whether this classifica- 
tion carries over in a meaningful way to the dynamics of 
spatially extended systems. 

We consider the set of Eqs. ([l])-([2| as canonical for two 
reasons. First, because the activator Eq. ([T]) has the sim- 
plest polynomial form of bistability. Note that Eq. ([T]) 
was for this reason originally suggested by Hodgkin and 
Huxley as the first mathematical model of the potas- 
sium dynamics in SD. It was published by Grafstein, who 
also provided experimental data supporting such a sim- 
ple reaction-diffusion scheme for the front dynamicP^. 
Second, the inhibitor Eq. ([2| has a linear rate function, 
in fact, the rate function is only of a function of the acti- 
vator u. This is the simplest inhibitor dynamics needed 
for pulse propagation. By neglecting an additional lin- 
ear term —^v in the inhibitor rate function, we limit 
the origin of excitability to the case of a supercritical 
Hopf bifurcation with subsequent canard explosion and 
avoid the bistable regime that exits in the subcritical 
case. The subcritical Hopf bifurcation occurs only in a 
narrow regime when 7 is close to 1 and /3 close to 0. We 
have tested some simulations with 7 = 0.5 with similar 
results. 

As a consequence of our assumptions about the model 
being in this canonical form, only two parameters exits, 
P which is associated with the threshold and the £, the 
time scale separation of activator and inhibitor dynam- 
ics. Of course, the choice of parameters can be quite 
different, a common choice is a in the cubic rate func- 
tion f{u) = u{u — a){u — a) but there are only two free 
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parameters or two equivalent groups of parameters. So 
there are the same bifurcations in the parameter planes 
(e, P) or (e, a), but to map the dynamics between equiva- 
lent groups of parameters might involve changes in time, 
space and concentrations scales. 

In particular the question of how the incidence of 
MA is reflected in the distance to the saddle-node bi- 
furcation, involves a measure on the parameter space, 
which we have suggested to get from pharmacokinetic- 
pharmacodynamic modelP^. 



Application to migraine pathophysiology 

The cause of the neurological symptoms in migraine 
with aura (MA) is the phenomenon of cortical spread- 
ing depression (SD)^^^^^^^. Whether SD is also a key 
to the subsequent headache phase is an open question, 
in particular, in cases of migraine without aura (MO ). If 
SD occurs in MO, it must remain clinically silenlP^or — 
by definition of diagnostic criteria — neurological symp- 
toms must last less than 5min. Of course, the transient 
nature of SD poses challenging problems in clinical ob- 
servability, in particular for objective measures by means 
of non-invasive imaging when clinical symptoms do not 
even indicate the aura phase with SD. The aura is usually, 
though not always, before the headache phase. Attacks 
observed with non-invasive imaging are usually triggered, 
which also could cause a trigger-specific bias. One well- 
documented case of a spontaneous migraine headache 
supports the contested notion of 'silent aura', because 
blood-fiow changes where observed that were most likely 
the result of SD^^ 

We suggest a qualitative congruence between the 
prevalence of MO and MA with the statistical properties 
we found in the transient response properties. We do not 
suggest that all MO attacks are related to SD nor that 
pain formation in MA is exclusively caused by SD. Rather 
that SD is one pathway of pain formation in MO and MA. 
We refer to this pathway as the "spreading depression" - 
theory of migraine^^. The "migraine generator" -theory 
(MG), a dysfunction in a central pattern generator in the 
brainstem that modulates the perception of pain, is for 
various reasons not less plausibl^. Some of the seem- 
ingly confiicting and controversially discussed evidence 
is probably resolved when one considered the basis of the 
classification of migraine. We currently have a symptom- 
based classification for migraine with possibly overlap- 
ping etiologies for individual subforms. In the light of 
an etiology-based classification with possibly overlapping 
symptoms the confiicts seem less puzzling to us. We also 
need to investigate an interplay of SD and MG, namely 
to which degree MG modulates pain traffic from SD gen- 
erated in the intracranial tissues. 

The cortex is not pain sensitive. There are detailed in- 
vestigations how SD in the cortex can cause pain via pain 
sensitive intracranial tissues and subsequent activation in 
the trigeminal nucleus caudalis in the brainstem^^ but 



dura 



dural sinuses 



sensory innervation 




small MIA 



cortex 



FIG. 9. Schematic representation of cross section of 
cortex, meninges and skull. The leptomeninges refer to 
the pia mater and arachnoid membrane. SD releases nox- 
ious substances with increased blood flow thought to diffuse 
outward. Activation of pain pathways can depend on MIA. 



^£| 53 | 54 [ rpj^^ qualitative congruence between the preva- 
lence of MO and MA with the statistical properties we 
found in the transient response properties is based on 
the following assumption on the geometrical layout: In 
the initial phase of cortical SD, with increased blood flow 
(hyperemic phase), a local release of noxious substances 
(ATP, glutamate, K+, H+) are thought to diffuse out- 
ward in the direction perpendicular to the cortex into 
"the leptomeninges resulting in activation of pial nocicep- 
tors, local neurogenic inflammation and the persistent ac- 
tivation of dural nociceptors which triggers the migraine 
headache"!^, but for issues concerning the blood brain 
barrier system cf.l^. If diffusion vertical to the affected 
area is critical, size and shape of this area should play 
a critical role, see Fig. [9j This suggests that SD waves 
activate nociceptive mechanisms dependent upon a suf- 
ficiently large instantaneously affected cortical area, i.e., 
large MIA. 

The aura phase on the other hand must clearly cor- 
relate with long and large enough cortical tissue being 
affected to notice the neurological deficits. In particular, 
because the very noticeable visual symptoms often start 
where the cortical magnification factor is large, so that 
only if they move into regions of lower magnification they 
get magnified by the reversed topographic mappin^^. 
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The seemingly contested notion of MO (migraine with- 
out aura) with silent aura is also resolved. 

A view at the connection between MIA and TAA as 
well as MIA and ED in our model is shown in Fig. 10 It 
shows that in the range of high MIA the average values 
for TAA and ED are becoming smaller. From Fig. [To] 
we can also read off that the range with the most events 
is in the regime of relatively high MIA (around 30) and 
significantly after the peak of ED resp. TAA. Moreover 
in the range with most events, the correlation coefficient 
r(MIA, ED) is always negative and the correlation coef- 
ficient r(MIA,TAA) is mostly negative. All these effects 
are stronger, the closer the control line is located to the 
saddle-node bifurcation. 

From these statistical correlations between MIA and 
TAA resp. ED and the distribution of the number of 
events, one could speculate that cases of MA are more 
rare and the quality of the headache in these cases might 
be less severe. This is exactly what has been reported in 
the medical literature^''. 

While the number of events with high ED and high 
TAA is inffuenced by the distance of the control line to 
the saddle node bifurcation, the number of events with 
high MIA is much lesser affected. So in a way, the dis- 
tance to the saddle node bifurcation controls the preva- 
lence of MA in our model, while the prevalence of MO is 
not much affected. 
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FIG. 10. Statistical analysis of output data. For all 

four pictures we took all data points with MIA in the interval 
[MIAiow, MIAiow+10] ("sliding window") and analyzed the 
connection with TAA (left column) and ED (right col- 
umn). In the (upper row), the average value is plotted 
with solid lines, the area of one standard deviation around 
this value is shaded. In the (lower row), the correlation co- 
efficient between MIA and the respective quantity is plotted. 
In all plots, the dotted lines indicates the number of events 
for the respective interval with MIAiow 



Inhibitory feedback and neurovascular coupling 

This naturally rises the question of the physiological 
origin of the inhibitory feedback control. The hyperemic 
phase engulfs large regions of the human cortex^^, while 
we suggest in this study the homeostatic breakdown di- 
rectly due to SD is much more limited in extent. A fast 
spreading increased neural activation in adjacent cortical 
areas could represent synaptic activation through feed- 
forward and feedback circuitry. This was suggested by 
Wilkinsori^. This would in turn extend the area of the 
hyperemic phase towards tissue that is not yet recruited 
into the SD state and this mechanism has therefore a 
neuroprotective effect by an increased blood flow which 
we mimic by the inhibitory mean field feedback. 

The coupling between neural activity and subsequent 
changes in cerebral blood flow, called neurovascular cou- 
pling, has a significant time delay in the order of seconds, 
which we ignore for the sake of simplicity in our model. 

Model-based control by neuromodulation 

We briefly discuss model-based control and means by 
which neuromodulation techniques may affect pathways 
of pain formation and the aura phase. 

The emerging transient patterns and their classifica- 
tion according to size and duration offer a model-based 
analysis of phase-dependent stimulation protocols for 



non-invasive neuromodulation devices, e.g. utilizing tran- 
scranial magnetic stimulation (TMS)^^, to intelligently 
target migraine. For instance, noise is a very effective 
method to drive the system back into the homogeneous 
steady state more quickly. In general, responses of non- 
linear systems to noise applied when the system is just 
before or past a saddle-node bifuraction are well stud- 
ied. Before the saddle-node on limit cycle bifuraction, 
the phenomenon of coherence resonance (CR) describes 
that a certain amount of noise makes responses most 
coherent^^. Behind the saddle-node bifurcation on a limit 
cycle the time the flows spend in the bottleneck region of 
the ghost is shortenecF^. However, noise would, accord- 
ing to our model, mainly positively affect ED and TAA, 
that is, the aura, while it could even worsen the headache, 
if applied early during the nucleation and growth process. 
Therefore, TMS using noise stimulation protocols, which 
is currently considered, should be applied only some time 
after first noticing aura symptoms. 

Headaches are not generally considered appropriate for 
invasive neurosurgical therapy, but when all else fails — 
preventives, abortives, and pain management — invasive 
brain stimulation techniques are also c onsidered, e.g. oc- 
cipital nerve stimulation (ONS)P^. So model-based 
control will become increasingly important. Also the im- 
portance of modeling related epileptic seizure dynamics 
as spatio-temporal transient patterns, has been suggested 
in a recent papei!^. Model-based control of Parkinson's 
disease, is already considered, yet Schiff remarks quite 
correctl^EH "It seems incredible that the tremendous 
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body of skill and knowledge of model-based control en- 
gineering has had so little impact on modern medicine. 
The timing is now propitious to propose fusing control 
theory with neural stimulation for the treatment of dy- 
namical brain disease." 

We suggest to consider migraine as a dynamical disease 
that could benefit from model-based control therapies. 



METHODS 

As a generic model for our excitable medium we use the 
well known FitzHugh-Nagumo equations^^, augmented 
by a diffusion term for the activator variable: 



du 

dv 
di 



= u 



= u 



1 



(1) 

(2) 



The parameter e separates the timescales of the dynamics 
of the activator u and the inhibitor v. £ is taken to be 
small. In the present work, we use a value of e = 0.04. 
The parameter /3 is a threshold value which determines 
from which activator level on the inhibitor concentration 
is rising. The local dynamics of ([2| (i.e. without the 
diffusion term) is oscillatory for \P\ < 1 and excitable 
for > 1. At = 1 the local dynamics undergo 
a supercritical Hopf-bifurcation. We choose a value of 
/3 = 1.1 throughout this work. To integrate ([2|, we used 
a simulation based on spectral methods^^ and adaptive 
timestepping. 

We define the (instantaneous) wave size as the area 
with activator level u over a certain threshold uq\ 



Sit) :-- 



11^ 



e {u{x, y, t) - '^threshold) dxdy , (3) 



where O is the Heaviside function and ^/threshold = 0- 

Equations ([2| are a paradigmatic model of an excitable 
medium^. It possesses a stable homogeneous solution as 
well as st able excited states (pulses, spirals or double 
spirals) cf !^^*^^* The boundary separating the basins of 
attraction of these types of solution consists of unstable 
so-called 'nucleation-solutions', which are areas of exci- 
tation which are traveling at uniform speed. The size in 
the sense of ([3| of these solutions is plotted against the 
parameter {3 in Fig. [3j To measure this line 9r, called 
the 'rotor boundary', we used a pseudo-continuation pro- 
cedure. Which is described below. 

Making the parameter j3 dependent on the wave size S 
adds a mean field control to the system. 



(4) 



where and /3o are control parameters. 

If the control line defined by Q intersects 9r, the point 
of intersection with higher S is stabilized, cf.^^. 

The aim of the present work is to shed light on the 
transient behavior, occurring when the control line Q is 
close to 9r but does not intersect it. 



To account for the imprecision in the rotor boundary 
of the simulation and the exact rotor boundary 9r, we 
measured the rotor boundary in our simulation using 
a pseudo-continuation procedure. For this, we set the 
control such that it intersects 9r and thus stabilize an 
otherwise unstable solution on it. Letting the simula- 
tion run until the system has stopped fluctuating, sav- 
ing the (/3,5')-pair, changing the control slightly and do- 
ing things over yields points of 9r in our system. From 
this measured 9r and the propagation boundary, inferred 
from continuation in ID, we chose 3 suitable control lines 
which were used for simulations in this work: 



K = 0.003 

(3o e [1.32,1.33,1.34] 



(5) 



For the purpose of visualization of the saddle-node bi- 
furcation occurring in the bifurcation diagram of the sys- 
tem with mean field control (right of Fig. |3j we fitted the 
measured branch of nucleation solutions to a function of 



the form f3 = a 



which also allows us to obtain 



cS+S'^ 

an approximation for the P value of the rotor boundary 

SRoo- 

As mentioned in section [j we need an appropriate sam- 
pling of initial conditions for ([2|, ideally being equidis- 
tantly sampled in some distribution. As was also men- 
tioned, the set of all initial conditions for this system 
does not — to our knowledge — carry a helpful mathemat- 
ical structure which allows us to achieve this aim easily. 
In order to attack this problem, we turned to the physi- 
ological origin for that we chose this model. 

A set of initial conditions should naturally reflect plau- 
sible spatial perturbations of the homogeneous steady 
state of the cortex. This can be achieved by defining lo- 
calized but spatially structured activity states on large- 
scales, i.e., on the order of millimeters. Such pattern 
are obtained from cortical features maps (see Fig. 11) 
by sampling three parameters {scaling^ depths and size) 
that define patches of lateral coupling in theses maps. A 
fourth parameter {excess) determines the amplitude of 
the perturbation. In the following, we first describe the 
rational behind using a cortical feature map and then the 
sampling. 

We focus on a cortical feature map in the primary vi- 
sual cortex (VI) called pinwheel map. VI is located at 
the occipital pole of the cerebral cortex and is the first 
region to process visual information from the eyes. Mi- 
graine aura symptoms often start there or nearby where 
similar feature maps exist. 

In VI, neurons within vertical columns (through the 
cortical layers) represent by their activity pattern edges, 
enlongated contours, and whole textures "seen" in the 
visual field. This representation has a distinct periodi- 
cally microstructured pattern: the pinwheel map. Neu- 
rons preferentially fire for edges with a given orientation 
and the preference changes continuously as a function 
of cortical location, except at singularities, the pinwheel 
centers, where the all the different orientations meetP^'^. 
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FIG. 11. To construct initial conditions from artificially gen- 
erated pinwheel maps, we first took such a pinwheel map with 
a certain scaling (upper left), then we chose a selection of ex- 
cited orientations by means of a gaussian. The width of the 
gaussian gives the selection depth (upper right). After that, 
we masked the result spatially with another gaussian distri- 
bution that is radially symmetric (lower right). The width of 
this gaussian gives the third parameter, the size of the pat- 
tern. Finally the result is scaled, giving rise to the fourth 
parameter, we called the excess and added to the activator 
variable in the homogeneous state (lower left). The inhibitor 
variable is put into the homogeneous state. 



Iso-orientation domains form continuous bands or 
patches around pinwheels and, on average, a region of 
about 1 mm^ (hypercolumn) will contain all possible ori- 
entation preferences. This topographical arrangement al- 
lows one hypercolumn to analyze all orientations com- 
ing from a small area in the visual field, but, as a con- 
sequence, the cortical representation of continuous con- 
tours in the visual field would be depicted in a patchy, 
discontinuous fashion''^. In general, spatially separated 
elements are bound together by short- and long-range lat- 
eral connections. While the strength of the local short- 
range connection within one hypercolum is a graded func- 
tion of cortical distance, mostly independent of relative 
orientation^^, long-range connections over several hyper- 
colums connect only iso-o rientation domains of similar 
orientation preferenc d^^ l ^^ l Even nearby regions, which 
are directly excitatory connected, have an inhibitory 
component through local inhibitory interneurons and this 
is likely be used to analyze angular visual features such 
as corners or T junction^^. 

Given the arguments above, we can now obtain local- 
ized yet spatially structured activity states on the scale 
we aim for as initial conditions by using iso-orientation 
domains that form continuous patches around pinwheels 
and extend in a discontinuous fashion over larger areas. 

IeP^ the authors analyzed the design principles that 



lie behind the columnar organization of the visual cor- 
tex. The precise design principles of this cortical organi- 
zation is governed by an annulus-like spectral structure 
in Fourier domain^^ , which is determined by mainly 
one parameter (scaling)^ that is, the annulus width. The 
parameter depth reflects the tuning properties of orien- 
tation preference or we can also interpret this as the 
range of orientation angles that we consider within the 
iso-orientation domain. The third parameter reflects the 
distance long-range coupling ranges before it significantly 
attenuates. 

These design principles can be exploited and a proce- 
dure can be designed to construct maps with the same 
properties. The constructed maps come very close to the 
maps found in brains of macaque monkeys (see^^ and 
references therein). 

To construct initial conditions from these maps we 
used a procedure that uses four control parameters and 
is visualized in Fig. 11 The details are as follows: A 



pinwheel map is a function that maps our twodimen- 
sional plane to the interval (— 7r/2, 7r/2]. We construct 
such a map using the procedure iiP^. During construc- 
tion, we can choose the scaling of the map. This is our 
first parameter. After constructing this map, by means 
of a gaussian, we choose a range of orientations that is 
excited. Mathematically speaking this is the concatena- 
tion of the gaussian distribution with the pinwheel map. 
This gives the next parameter, namely the width of the 
gaussian that selects the angles, we call that parameter 
the depth. The next step is to constrain the generated 
pattern spatially by multiplication with another gaussian 
which is defined on the plane P and chosen to be rota- 
tionally symmetric. The width of this gaussian gives rise 
to the third parameter, the size of the pattern. Finally 
we multiply the pattern by a certain amplitude, which 
is chosen such that the integral of the pattern over the 
plane gives a chosen number, which constitutes the fourth 
parameter, we called the excess. 

Finally, initial conditions are generated by setting the 
plane to the (stable) homogeneous state and then adding 
the generated pattern to the activator variable u. 

In a first run, we scanned the space spanned by the four 
parameters coarsely. We used the marginal distributions 
of the number of solutions with ED> with respect to 
the parameters to decide how densely to sample the pa- 
rameter space in the final run. 
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